Assignment 7

Setup

knitr::opts_chunk$set(echo = TRUE, comment = "#>", dpi = 300)

for (f in list.files(here::here("src"), pattern = "R$", full.names = TRUE)) {
  source(f)
}

library(rstan)
library(tidybayes)
library(magrittr)
library(tidyverse)

theme_set(theme_classic() + theme(strip.background = element_blank()))

options(mc.cores = 2)
rstan_options(auto_write = TRUE)

drowning <- aaltobda::drowning
factory <- aaltobda::factory

Assignment 7

1. Linear model: drowning data with Stan

The provided data drowning in the ‘aaltobda’ package contains the number of people who died from drowning each year in Finland 1980–2019. A statistician is going to fit a linear model with Gaussian residual model to these data using time as the predictor and number of drownings as the target variable. She has two objective questions:

  1. What is the trend of the number of people drowning per year? (We would plot the histogram of the slope of the linear model.)
  2. What is the prediction for the year 2020? (We would plot the histogram of the posterior predictive distribution for the number of people drowning at \(\tilde{x} = 2020\).)

Corresponding Stan code is provided in Listing 1. However, it is not entirely correct for the problem. First, there are three mistakes. Second, there are no priors defined for the parameters. In Stan, this corresponds to using uniform priors.

a) Find the three mistakes in the code and fix them. Report the original mistakes and your fixes clearly in your report. Include the full corrected Stan code in your report.

  1. Declaration of sigma on line 10 should be real<lower=0>.
  2. Missing semicolon at the end of line 16.
  3. On line 19, the prediction on new data does not use the new data in xpred. This has been changed to real ypred = normal_rng(alpha + beta*xpred, sigma);.

Below is a copy of the final model. The full Stan file is at models/assignment07-drownings.stan.

data {
  int<lower=0> N;  // number of data points
  vector[N] x;     // observation year
  vector[N] y;     // observation number of drowned
  real xpred;      // prediction year
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;  // fix: 'upper' should be 'lower'
}
transformed parameters {
  vector[N] mu = alpha + beta*x;
}
model {
  y ~ normal(mu, sigma);  // fix: missing semicolor
}
generated quantities {
  real ypred = normal_rng(alpha + beta*xpred, sigma);  // fix: use `xpred`
}

b) Determine a suitable weakly-informative prior \(\text{Normal}(0,\sigma_\beta)\) for the slope \(\beta\). It is very unlikely that the mean number of drownings changes more than 50 % in one year. The approximate historical mean yearly number of drownings is 138. Hence, set \(\sigma_\beta\) so that the following holds for the prior probability for \(\beta\): \(Pr(−69 < \beta < 69) = 0.99\). Determine suitable value for \(\sigma_\beta\) and report the approximate numerical value for it.

x <- rnorm(1e5, 0, 26)
print(mean(-69 < x & x < 69))

#> [1] 0.99179
plot_single_hist(x, alpha = 0.5, color = "black") + geom_vline(xintercept = c(-69, 69)) + labs(x = "beta")

c) Using the obtained σβ, add the desired prior in the Stan code.

From some trial and error, it seems that a prior of \(\text{Normal}(0, 26)\) should work. I have added this prior distribution to beta in the model at line 17.

beta ~ normal(0, 26);   // prior on `beta`

d) In a similar way, add a weakly informative prior for the intercept alpha and explain how you chose the prior.

To use the year directly as the values for \(x\) would lead to a massive value of \(\alpha\) because the values for \(x\) range from 1980 to 2019. Thus, it would be advisable to first center the year, meaning at the prior distribution for \(\alpha\) can be centered around the average of the number of drownings per year and a standard deviation near that of the actual number of drownings.

head(drowning)

#>   year drownings
#> 1 1980       149
#> 2 1981       127
#> 3 1982       139
#> 4 1983       141
#> 5 1984       122
#> 6 1985       120
print(mean(drowning$drownings))

#> [1] 134.35
print(sd(drowning$drownings))

#> [1] 28.48441

Therefore, I add the prior \(\text{Normal}(135, 50)\) to \(\alpha\) on line 16.

alpha ~ normal(135, 50);   // prior on `alpha`
data <- list(
  N = nrow(drowning),
  x = drowning$year - mean(drowning$year),
  y = drowning$drownings,
  xpred = 2020 - mean(drowning$year)
)
drowning_model <- stan(
  here::here("models", "assignment07-drownings.stan"),
  data = data
)

variable_post <- spread_draws(drowning_model, alpha, beta) %>%
  pivot_longer(c(alpha, beta), names_to = "variable", values_to = "value")
head(variable_post)

#> # A tibble: 6 × 5
#>   .chain .iteration .draw variable   value
#>    <int>      <int> <int> <chr>      <dbl>
#> 1      1          1     1 alpha    133.
#> 2      1          1     1 beta      -0.863
#> 3      1          2     2 alpha    134.
#> 4      1          2     2 beta      -1.14
#> 5      1          3     3 alpha    133.
#> 6      1          3     3 beta      -1.28
variable_post %>%
  ggplot(aes(x = .iteration, y = value, color = factor(.chain))) +
  facet_grid(rows = vars(variable), scales = "free_y") +
  geom_path(alpha = 0.5) +
  scale_x_continuous(expand = expansion(c(0, 0))) +
  scale_y_continuous(expand = expansion(c(0.02, 0.02))) +
  labs(x = "iteration", y = "value", color = "chain")


variable_post %>%
  ggplot(aes(x = value)) +
  facet_grid(cols = vars(variable), scales = "free_x") +
  geom_histogram(color = "black", alpha = 0.3, bins = 30) +
  scale_x_continuous(expand = expansion(c(0.02, 0.02))) +
  scale_y_continuous(expand = expansion(c(0, 0.02)))


spread_draws(drowning_model, ypred) %$%
  plot_single_hist(ypred, alpha = 0.3, color = "black") +
  labs(x = "predicted number of drownings in 2020")


red <- "#C34E51"

bayestestR::describe_posterior(drowning_model, ci = 0.89, test = c()) %>%
  as_tibble() %>%
  filter(str_detect(Parameter, "mu")) %>%
  select(Parameter, Median, CI_low, CI_high) %>%
  janitor::clean_names() %>%
  mutate(idx = row_number()) %>%
  left_join(drowning %>% mutate(idx = row_number()), by = "idx") %>%
  ggplot(aes(x = year)) +
  geom_point(aes(y = drownings), data = drowning, color = "#4C71B0") +
  geom_line(aes(y = median), color = red, size = 1.2) +
  geom_smooth(
    aes(y = ci_low),
    method = "loess",
    formula = "y ~ x",
    linetype = 2,
    se = FALSE,
    color = red,
    size = 1
  ) +
  geom_smooth(
    aes(y = ci_high),
    method = "loess",
    formula = "y ~ x",
    linetype = 2,
    se = FALSE,
    color = red,
    size = 1
  ) +
  labs(x = "year", y = "number of drownings (mean ± 89% CI)")

2. Hierarchical model: factory data with Stan

The factory data in the ‘aaltobda’ package contains quality control measurements from 6 machines in a factory (units of the measurements are irrelevant here). In the data file, each column contains the measurements for a single machine. Quality control measurements are expensive and time-consuming, so only 5 measurements were done for each machine. In addition to the existing machines, we are interested in the quality of another machine (the seventh machine).

For this problem, you’ll use the following Gaussian models:

As in the model described in the book, use the same measurement standard deviation \(\sigma\) for all the groups in the hierarchical model. In the separate model, however, use separate measurement standard deviation \(\sigma_j\) for each group \(j\). You should use weakly informative priors for all your models.

Complete the following questions for each of the three models (separate, pooled, hierarchical).

a) Describe the model with mathematical notation. Also describe in words the difference between the three models.

Separate model: The separate model is described below where each machine has its own centrality \(\mu\) and dispersion \(\sigma\) parameters that do not influence the parameters of the other machines.

\[ y_{ij} \sim N(\mu_j, \sigma_j) \\ \mu_j \sim N(0, 1) \\ \sigma_j \sim \text{Inv-}\chi^2(10) \]

Pooled model: The pool model is described below where there is no distinction between the models but instead a single set of parameters for all of the data.

\[ y_{i} \sim N(\mu, \sigma) \\ \mu \sim N(0, 1) \\ \sigma \sim \text{Inv-}\chi^2(10) \]

Hierarchical model: The hierarchical model is described below where each machine has its own centrality \(\mu\) parameter which are linked through a hyper-prior distribution from which they are drawn. The machines will all share a common dispersion paramete \(\sigma\)

\[ y_{ij} \sim N(\mu_j, \sigma_j) \\ \mu_j \sim N(\alpha, \tau) \\ \alpha \sim N(0, 1) \\ \tau \sim \text{HalfNormal}(2.5) \\ \sigma \sim \text{Inv-}\chi^2(10) \]

The separate model is effectively building a different linear model for each machine where as the pooled model treats all the measurements as coming from the same model. The hierarchical model is treating the machines as having come from a single, shared distribution.

b) Implement the model in Stan and include the code in the report. Use weakly informative priors for all your models.

print_model_code <- function(path) {
  for (l in readLines(path)) {
    cat(l, "\n")
  }
}

Separate model

separate_model_code <- here::here(
  "models", "assignment07_factories_separate.stan"
)
print_model_code(separate_model_code)

#> data {
#>   int<lower=0> N;  // number of data points per machine
#>   int<lower=0> J;  // number of machines
#>   vector[J] y[N];  // quality control data points
#> }
#>
#> parameters {
#>   vector[J] mu;
#>   vector<lower=0>[J] sigma;
#> }
#>
#> model {
#>   // priors
#>   for (j in 1:J) {
#>     mu[j] ~ normal(100, 10);
#>     sigma[j] ~ inv_chi_square(5);
#>   }
#>
#>   // likelihood
#>   for (j in 1:J){
#>     y[,j] ~ normal(mu[j], sigma[j]);
#>   }
#> }
#>
#> generated quantities {
#>   // Compute the predictive distribution for the sixth machine.
#>   real y6pred;
#>   vector[J] log_lik[N];
#>
#>   y6pred = normal_rng(mu[6], sigma[6]);
#>
#>   for (j in 1:J) {
#>     for (n in 1:N) {
#>       log_lik[n,j] = normal_lpdf(y[n,j] | mu[j], sigma[j]);
#>     }
#>   }
#> }
separate_model_data <- list(
  y = factory,
  N = nrow(factory),
  J = ncol(factory)
)
separate_model <- rstan::stan(
  separate_model_code,
  data = separate_model_data,
  verbose = FALSE,
  refresh = 0
)
knitr::kable(
  bayestestR::describe_posterior(separate_model, ci = 0.89, test = NULL),
  digits = 3
)
Parameter Median CI CI_low CI_high ESS Rhat
31 mu[1] 85.109 0.89 73.881 97.273 3909.502 1.001
32 mu[2] 105.207 0.89 97.768 111.954 3279.706 1.001
33 mu[3] 90.084 0.89 83.087 97.754 3561.168 1.000
34 mu[4] 110.779 0.89 105.953 115.134 3882.334 1.000
35 mu[5] 91.402 0.89 85.182 97.764 3573.824 0.999
36 mu[6] 90.748 0.89 81.428 100.556 4594.641 1.000
37 y6pred 90.502 0.89 60.694 119.899 3896.319 0.999
1 log_lik[1,1] -3.866 0.89 -4.431 -3.379 2176.770 1.004
7 log_lik[2,1] -3.985 0.89 -4.452 -3.537 5132.386 1.000
13 log_lik[3,1] -3.985 0.89 -4.452 -3.537 4647.762 1.000
19 log_lik[4,1] -6.257 0.89 -8.205 -4.776 4641.331 1.000
25 log_lik[5,1] -4.392 0.89 -5.077 -3.762 5477.005 1.000
2 log_lik[1,2] -4.020 0.89 -4.880 -3.252 6389.885 0.999
8 log_lik[2,2] -3.346 0.89 -3.911 -2.869 2972.969 1.002
14 log_lik[3,2] -3.684 0.89 -4.334 -3.040 2158.970 1.000
20 log_lik[4,2] -3.290 0.89 -3.768 -2.792 2411.683 1.001
26 log_lik[5,2] -4.991 0.89 -6.987 -3.639 6650.017 0.999
3 log_lik[1,3] -3.907 0.89 -4.697 -3.284 4032.901 0.999
9 log_lik[2,3] -3.419 0.89 -3.893 -2.991 3290.788 1.000
15 log_lik[3,3] -3.388 0.89 -3.884 -2.976 2972.969 1.002
21 log_lik[4,3] -3.425 0.89 -3.985 -2.933 3733.757 1.000
27 log_lik[5,3] -5.682 0.89 -7.642 -4.110 2285.753 1.000
4 log_lik[1,4] -3.290 0.89 -3.996 -2.682 4037.383 1.000
10 log_lik[2,4] -3.695 0.89 -4.710 -2.892 5346.134 0.999
16 log_lik[3,4] -3.206 0.89 -3.872 -2.589 4088.675 1.000
22 log_lik[4,4] -3.786 0.89 -5.056 -2.872 7472.336 1.000
28 log_lik[5,4] -3.206 0.89 -3.872 -2.589 2017.422 1.000
5 log_lik[1,5] -4.143 0.89 -5.153 -3.231 2362.561 1.000
11 log_lik[2,5] -3.409 0.89 -3.953 -2.920 5353.519 1.000
17 log_lik[3,5] -4.029 0.89 -5.114 -3.229 5477.005 1.000
23 log_lik[4,5] -4.143 0.89 -5.153 -3.231 4032.005 1.000
29 log_lik[5,5] -3.189 0.89 -3.687 -2.697 3772.266 1.000
6 log_lik[1,6] -5.914 0.89 -7.864 -4.593 4952.273 1.000
12 log_lik[2,6] -3.769 0.89 -4.235 -3.311 6883.850 1.000
18 log_lik[3,6] -4.145 0.89 -4.692 -3.612 4037.383 1.000
24 log_lik[4,6] -4.134 0.89 -4.721 -3.579 2736.658 1.000
30 log_lik[5,6] -3.965 0.89 -4.418 -3.511 3743.014 1.000

Pooled model

pooled_model_code <- here::here("models", "assignment07_factories_pooled.stan")
print_model_code(pooled_model_code)

#> data {
#>   int<lower=0> N;  // number of data points
#>   vector[N] y;     // machine quality control data
#> }
#>
#> parameters {
#>   real mu;
#>   real<lower=0> sigma;
#> }
#>
#> model {
#>   // priors
#>   mu ~ normal(100, 10);
#>   sigma ~ inv_chi_square(5);
#>
#>   // likelihood
#>   y ~ normal(mu, sigma);
#> }
#>
#> generated quantities {
#>   real ypred;
#>   vector[N] log_lik;
#>
#>   ypred = normal_rng(mu, sigma);
#>
#>   for (i in 1:N)
#>     log_lik[i] = normal_lpdf(y[i] | mu, sigma);
#>
#> }
pooled_model_data <- list(
  y = unname(unlist(factory)),
  N = length(unlist(factory))
)
pooled_model <- rstan::stan(
  pooled_model_code,
  data = pooled_model_data,
  verbose = FALSE,
  refresh = 0
)

knitr::kable(
  bayestestR::describe_posterior(pooled_model, ci = 0.89, test = NULL),
  digits = 3
)
Parameter Median CI CI_low CI_high ESS Rhat
31 mu 93.685 0.89 88.720 98.481 2272.795 1.004
32 ypred 93.585 0.89 63.431 120.626 3581.065 1.002
1 log_lik[1] -3.980 0.89 -4.194 -3.752 2282.741 1.002
12 log_lik[2] -3.794 0.89 -4.012 -3.594 2198.446 1.001
23 log_lik[3] -3.794 0.89 -4.012 -3.594 2198.446 1.001
25 log_lik[4] -7.526 0.89 -8.930 -6.015 3059.601 1.001
26 log_lik[5] -4.956 0.89 -5.429 -4.472 3401.936 1.002
27 log_lik[6] -4.686 0.89 -5.143 -4.338 2450.418 1.003
28 log_lik[7] -4.183 0.89 -4.451 -3.954 3011.015 1.003
29 log_lik[8] -4.471 0.89 -4.825 -4.159 2470.049 1.003
30 log_lik[9] -3.974 0.89 -4.193 -3.773 2975.150 1.002
2 log_lik[10] -3.863 0.89 -4.082 -3.657 2160.166 1.001
3 log_lik[11] -3.884 0.89 -4.103 -3.693 2753.815 1.001
4 log_lik[12] -3.791 0.89 -4.006 -3.587 2228.389 1.000
5 log_lik[13] -3.794 0.89 -4.012 -3.594 2198.446 1.001
6 log_lik[14] -3.888 0.89 -4.098 -3.672 2178.957 1.002
7 log_lik[15] -4.956 0.89 -5.429 -4.472 3401.936 1.002
8 log_lik[16] -4.009 0.89 -4.242 -3.812 2981.695 1.002
9 log_lik[17] -4.851 0.89 -5.338 -4.416 2477.500 1.002
10 log_lik[18] -4.611 0.89 -4.999 -4.239 2437.549 1.003
11 log_lik[19] -3.912 0.89 -4.127 -3.715 2854.113 1.001
13 log_lik[20] -4.611 0.89 -4.999 -4.239 2437.549 1.003
14 log_lik[21] -4.146 0.89 -4.400 -3.912 2410.994 1.003
15 log_lik[22] -3.809 0.89 -4.017 -3.603 2422.870 1.000
16 log_lik[23] -3.941 0.89 -4.156 -3.739 2952.412 1.002
17 log_lik[24] -4.146 0.89 -4.400 -3.912 2410.994 1.003
18 log_lik[25] -3.794 0.89 -4.012 -3.594 2198.446 1.001
19 log_lik[26] -5.986 0.89 -6.861 -5.146 3232.831 1.001
20 log_lik[27] -3.794 0.89 -4.012 -3.594 2198.446 1.001
21 log_lik[28] -3.974 0.89 -4.193 -3.773 2975.150 1.002
22 log_lik[29] -4.251 0.89 -4.518 -3.987 2627.063 1.003
24 log_lik[30] -3.859 0.89 -4.077 -3.668 2658.622 1.001

Hierarchical model

hierarchical_model_code <- here::here(
  "models", "assignment07_factories_hierarchical.stan"
)
print_model_code(hierarchical_model_code)

#> data {
#>   int<lower=0> N;  // number of data points per machine
#>   int<lower=0> J;  // number of machines
#>   vector[J] y[N];  // quality control data points
#> }
#>
#> parameters {
#>   vector[J] mu;
#>   real<lower=0> sigma;
#>   real alpha;
#>   real<lower=0> tau;
#> }
#>
#> model {
#>   // hyper-priors
#>   alpha ~ normal(100, 10);
#>   tau ~ normal(0, 10);
#>
#>   // priors
#>   mu ~ normal(alpha, tau);
#>   sigma ~ inv_chi_square(5);
#>
#>   // likelihood
#>   for (j in 1:J){
#>     y[,j] ~ normal(mu[j], sigma);
#>   }
#> }
#>
#> generated quantities {
#>   // Compute the predictive distribution for the sixth machine.
#>   real y6pred;  // Leave for compatibility with earlier assignments.
#>   vector[J] ypred;
#>   real mu7pred;
#>   real y7pred;
#>   vector[J] log_lik[N];
#>
#>   y6pred = normal_rng(mu[6], sigma);
#>   for (j in 1:J) {
#>     ypred[j] = normal_rng(mu[j], sigma);
#>   }
#>
#>   mu7pred = normal_rng(alpha, tau);
#>   y7pred = normal_rng(mu7pred, sigma);
#>
#>   for (j in 1:J) {
#>     for (n in 1:N) {
#>       log_lik[n,j] = normal_lpdf(y[n,j] | mu[j], sigma);
#>     }
#>   }
#> }
hierarchical_model_data <- list(
  y = factory,
  N = nrow(factory),
  J = ncol(factory)
)
hierarchical_model <- rstan::stan(
  hierarchical_model_code,
  data = hierarchical_model_data,
  verbose = FALSE,
  refresh = 0
)
knitr::kable(
  bayestestR::describe_posterior(hierarchical_model, ci = 0.89, test = NULL),
  digits = 3
)
Parameter Median CI CI_low CI_high ESS Rhat
32 mu[1] 81.501 0.89 71.585 91.564 1551.901 1.002
33 mu[2] 102.585 0.89 93.157 112.082 1282.943 1.001
34 mu[3] 89.815 0.89 80.625 98.622 3179.351 1.000
35 mu[4] 106.492 0.89 96.550 117.140 884.281 1.003
36 mu[5] 91.400 0.89 82.756 100.395 3969.285 0.999
37 mu[6] 88.474 0.89 80.427 97.934 3071.141 1.000
1 alpha 94.328 0.89 87.300 102.025 2965.380 1.000
39 tau 10.590 0.89 4.483 17.607 805.385 1.003
40 y6pred 87.953 0.89 63.999 112.354 4313.435 1.000
42 ypred[1] 81.718 0.89 58.608 108.899 3513.350 1.000
43 ypred[2] 102.759 0.89 79.005 128.502 3067.405 1.001
44 ypred[3] 89.773 0.89 67.272 115.928 4153.920 0.999
45 ypred[4] 106.136 0.89 79.471 130.586 3296.583 1.000
46 ypred[5] 91.376 0.89 65.793 115.846 3920.138 0.999
47 ypred[6] 88.446 0.89 62.515 111.900 3954.587 1.001
38 mu7pred 94.340 0.89 73.758 113.571 3499.233 1.001
41 y7pred 95.018 0.89 65.644 126.849 3759.822 1.000
2 log_lik[1,1] -3.656 0.89 -3.972 -3.396 1420.286 1.003
8 log_lik[2,1] -3.881 0.89 -4.474 -3.487 1675.051 1.000
14 log_lik[3,1] -3.881 0.89 -4.474 -3.487 2963.128 1.001
20 log_lik[4,1] -6.739 0.89 -8.435 -4.876 749.030 1.006
26 log_lik[5,1] -4.117 0.89 -4.826 -3.495 3809.548 0.999
3 log_lik[1,2] -4.106 0.89 -4.798 -3.542 4462.815 1.000
9 log_lik[2,2] -3.712 0.89 -4.128 -3.368 3327.369 1.000
15 log_lik[3,2] -3.914 0.89 -4.521 -3.474 895.217 1.002
21 log_lik[4,2] -3.642 0.89 -3.953 -3.375 1689.538 1.002
27 log_lik[5,2] -4.197 0.89 -5.021 -3.647 1023.154 1.004
4 log_lik[1,3] -3.917 0.89 -4.485 -3.478 2751.946 1.001
10 log_lik[2,3] -3.663 0.89 -3.961 -3.391 2328.978 1.002
16 log_lik[3,3] -3.653 0.89 -3.945 -3.395 3327.369 1.000
22 log_lik[4,3] -3.665 0.89 -4.000 -3.385 1293.150 1.001
28 log_lik[5,3] -4.869 0.89 -5.984 -3.930 1608.427 1.002
5 log_lik[1,4] -3.659 0.89 -3.967 -3.389 783.916 1.005
11 log_lik[2,4] -3.980 0.89 -4.688 -3.462 3975.663 0.999
17 log_lik[3,4] -3.819 0.89 -4.424 -3.399 3237.502 1.000
23 log_lik[4,4] -3.696 0.89 -4.025 -3.390 4136.744 1.000
29 log_lik[5,4] -3.819 0.89 -4.424 -3.399 785.243 1.003
6 log_lik[1,5] -3.980 0.89 -4.569 -3.503 1682.013 1.000
12 log_lik[2,5] -3.709 0.89 -4.044 -3.405 1917.116 1.003
18 log_lik[3,5] -3.940 0.89 -4.490 -3.485 3809.548 0.999
24 log_lik[4,5] -3.980 0.89 -4.569 -3.503 3077.755 1.000
30 log_lik[5,5] -3.639 0.89 -3.904 -3.373 1843.842 1.001
7 log_lik[1,6] -6.033 0.89 -7.671 -4.724 2583.433 1.000
13 log_lik[2,6] -3.665 0.89 -3.935 -3.369 4406.506 1.000
19 log_lik[3,6] -4.192 0.89 -4.936 -3.608 783.916 1.005
25 log_lik[4,6] -3.922 0.89 -4.476 -3.475 1347.074 1.002
31 log_lik[5,6] -3.933 0.89 -4.476 -3.489 3391.934 1.001

c) Using the model (with weakly informative priors) report, comment on and, if applicable, plot histograms for the following distributions:

  1. the posterior distribution of the mean of the quality measurements of the sixth machine.
  2. the predictive distribution for another quality measurement of the sixth machine.
  3. the posterior distribution of the mean of the quality measurements of the seventh machine.
plot_hist_mean_of_sixth <- function(vals) {
  plot_single_hist(vals, bins = 30, color = "black", alpha = 0.3) +
    labs(x = "mean of 6th machine", y = "posterior density")
}

plot_hist_sixth_predictions <- function(vals) {
  plot_single_hist(vals, bins = 30, color = "black", alpha = 0.3) +
    labs(x = "posterior predictions for 6th machine", y = "posterior density")
}

plot_hist_mean_of_seventh <- function(vals) {
  plot_single_hist(vals, bins = 30, color = "black", alpha = 0.3) +
    labs(x = "mean of 7thth machine", y = "posterior density")
}

Separate model

plot_hist_mean_of_sixth(rstan::extract(separate_model)$mu[, 6])


plot_hist_sixth_predictions(rstan::extract(separate_model)$y6pred)

It is not possible to estimate the posterior for the mean of some new 7th machine because all machines are treated separately.

Pooled model

plot_hist_mean_of_sixth(rstan::extract(pooled_model)$mu)


plot_hist_sixth_predictions(rstan::extract(pooled_model)$ypred)

The predicted mean for a new machine is the same as the pooled mean \(mu\).

plot_hist_mean_of_seventh(rstan::extract(pooled_model)$mu)

Hierarchical model

plot_hist_mean_of_sixth(rstan::extract(hierarchical_model)$mu[, 6])


plot_hist_sixth_predictions(rstan::extract(hierarchical_model)$y6pred)


plot_hist_mean_of_seventh(rstan::extract(hierarchical_model)$mu7pred)

d) Report the posterior expectation for \(\mu_1\) with a 90% credible interval but using a \(\text{Normal}(0,10)\) prior for the \(\mu\) parameter(s) and a \(\text{Gamma}(1,1)\) prior for the \(\sigma\) parameter(s). For the hierarchical model, use the \(\text{Normal}(0, 10)\) and \(\text{Gamma}(1, 1)\) as hyper-priors.

(I’m going to skip this one, but come back to it if it is needed for future assignments.)

sessionInfo()

#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur 11.6
#>
#> Matrix products: default
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats     graphics  grDevices datasets  utils     methods
#> [7] base
#>
#> other attached packages:
#>  [1] magrittr_2.0.1       bayesplot_1.8.1      aaltobda_0.3.1
#>  [4] tidybayes_3.0.1      rstan_2.21.2         StanHeaders_2.21.0-7
#>  [7] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.7
#> [10] purrr_0.3.4          readr_2.0.1          tidyr_1.1.3
#> [13] tibble_3.1.3         ggplot2_3.3.5        tidyverse_1.3.1
#> [16] patchwork_1.1.1      ggtext_0.1.1         glue_1.4.2
#>
#> loaded via a namespace (and not attached):
#>   [1] colorspace_2.0-2     ggridges_0.5.3       ellipsis_0.3.2
#>   [4] rsconnect_0.8.24     rprojroot_2.0.2      snakecase_0.11.0
#>   [7] markdown_1.1         fs_1.5.0             gridtext_0.1.4
#>  [10] rstudioapi_0.13      farver_2.1.0         bit64_4.0.5
#>  [13] svUnit_1.0.6         fansi_0.5.0          mvtnorm_1.1-2
#>  [16] lubridate_1.7.10     xml2_1.3.2           codetools_0.2-18
#>  [19] splines_4.1.2        downlit_0.2.1        knitr_1.33
#>  [22] jsonlite_1.7.2       LaplacesDemon_16.1.6 broom_0.7.9
#>  [25] dbplyr_2.1.1         ggdist_3.0.0         compiler_4.1.2
#>  [28] httr_1.4.2           backports_1.2.1      assertthat_0.2.1
#>  [31] Matrix_1.3-4         cli_3.0.1            htmltools_0.5.1.1
#>  [34] prettyunits_1.1.1    tools_4.1.2          coda_0.19-4
#>  [37] gtable_0.3.0         posterior_1.1.0      V8_3.4.2
#>  [40] Rcpp_1.0.7           cellranger_1.1.0     jquerylib_0.1.4
#>  [43] vctrs_0.3.8          nlme_3.1-153         insight_0.14.4
#>  [46] tensorA_0.36.2       xfun_0.25            ps_1.6.0
#>  [49] rvest_1.0.1          lifecycle_1.0.0      renv_0.14.0
#>  [52] MASS_7.3-54          scales_1.1.1         vroom_1.5.4
#>  [55] clisymbols_1.2.0     hms_1.1.0            parallel_4.1.2
#>  [58] inline_0.3.19        RColorBrewer_1.1-2   yaml_2.2.1
#>  [61] curl_4.3.2           gridExtra_2.3        loo_2.4.1
#>  [64] sass_0.4.0           distill_1.2          stringi_1.7.3
#>  [67] bayestestR_0.11.0    highr_0.9            checkmate_2.0.0
#>  [70] pkgbuild_1.2.0       rlang_0.4.11         pkgconfig_2.0.3
#>  [73] matrixStats_0.61.0   distributional_0.2.2 evaluate_0.14
#>  [76] lattice_0.20-45      labeling_0.4.2       bit_4.0.4
#>  [79] processx_3.5.2       tidyselect_1.1.1     here_1.0.1
#>  [82] plyr_1.8.6           bookdown_0.23        R6_2.5.0
#>  [85] generics_0.1.0       DBI_1.1.1            pillar_1.6.2
#>  [88] haven_2.4.3          withr_2.4.2          mgcv_1.8-38
#>  [91] datawizard_0.2.1     abind_1.4-5          janitor_2.1.0
#>  [94] modelr_0.1.8         crayon_1.4.1         arrayhelpers_1.1-0
#>  [97] utf8_1.2.2           tzdb_0.1.2           rmarkdown_2.10
#> [100] isoband_0.2.5
#>  [ reached getOption("max.print") -- omitted 10 entries ]